%% Stock Return Extrapolation, Option Prices, and Variance Risk Premium
% By Adem Atmaz
% This file generates Table 6:Effects for the options: Model vs. evidence
% It also requires 2 function files: RN_OptionPrices.m and RN_OptionPsi.m
% NOTE: Typically it takes about 1 HOUR to run

clear all;  clc;
format compact
format

% Parameter Values follow from  "Atmaz_Table_4_ParameterValues.m"
alpha=0.50;
theta=0.25;
mu=0.0212;
kappa=0.3567;
zeta=0.0044;
sigma=0.1625;
rho=0;
r=0.0056;

D0=100;
a=2.75;
vec_theta=[0, 0.25, 0.50, 0.75];
N_th=length(vec_theta);


% Option contract quantities
Tt=2/12;       % Option Maturity is 2 months
K_OTM_P=0.96;  %  OTM put moneyness
K_OTM_C=1.04;  %  OTM call moneyness
K_ATM=1.00;    %  ATM put and call moneyness

N_sim=1000;     % Number of simulations
T=2/12;         % simulation horizon 2 months as in Amin
h = 1/252;      % Use 1/252 for daily returns
t = 0:h:T;      % Time vector from 0 to T
J=length(t);

% Initilization of matrices
Put_OTM_Pos_IV=NaN(N_sim,N_th);
Put_OTM_Neg_IV=NaN(N_sim,N_th);
Put_ATM_Pos_IV=NaN(N_sim,N_th);
Put_ATM_Neg_IV=NaN(N_sim,N_th);
Call_OTM_Pos_IV=NaN(N_sim,N_th);
Call_OTM_Neg_IV=NaN(N_sim,N_th);
Call_ATM_Pos_IV=NaN(N_sim,N_th);
Call_ATM_Neg_IV=NaN(N_sim,N_th);

vec_D = zeros(1,J); % Allocate output vector
vec_V = zeros(1,J); % Allocate output vector
vec_X = zeros(1,J); % Allocate output vector

tic
for th=1:N_th
    theta=vec_theta(th)   % Vary theta
    % Implied Quantities after fixing theta
    M_DELTA=((kappa+(rho*sigma)*(1+theta))^2)-((2+theta)*(1+theta)*sigma^2);
    M_c=theta/(alpha*(1+theta));
    M_b=-(2+theta)/(kappa+((rho*sigma)*(1+theta))+sqrt(M_DELTA));
    M_a=a;
    M_sigma_1=(1+(rho*sigma*M_b))/(1-alpha*M_c);
    M_sigma_2=(sqrt(1-rho^2)*sigma*M_b)/(1-alpha*M_c);
    M_Var_con=(M_sigma_1^2)+(M_sigma_2^2);
    M_Cov_con=(rho+(sigma*M_b))/(1-alpha*M_c);
    M_corr_sv=M_Cov_con/sqrt(M_Var_con);    
    M_zeta_v=zeta*M_Var_con;
    M_kappa_v=kappa+(sigma*M_Cov_con);
    M_eta_v=(theta/M_b)*(1-(alpha*M_c))*M_Var_con;
    M_sigma_v=sigma*sqrt(M_Var_con);
    M_rho_sv=M_Cov_con/sqrt(M_Var_con);
    % long-run mean (steady-state) values of processes
    D0=100;  % initial fundamental value
    V0=zeta/kappa; % long run mean of  fundamental  variance
    v0=M_zeta_v/kappa;  % LRM of stock return variance v
    X0=(r+(0.5*v0))/(1+theta); % LRM of sentiment X
    % Initalization of each path
    vec_D(1)=D0;
    vec_V(1)=V0;
    vec_X(1)=X0;    
    
    % ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    %  ####################### Option Prices #################################
    %  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
    
    rng('default');
    Z1_Mat=randn(N_sim,J);
    Z2_Mat=randn(N_sim,J);    
    for path=1:N_sim
        path        
        Z1=Z1_Mat(path,:);
        Z2=Z2_Mat(path,:);
        for j = 1:J-1
            vec_V(j+1) =  max(0, vec_V(j))+((zeta-(kappa* max(0, vec_V(j))))*h)...
                +(sigma*rho*sqrt( max(0, vec_V(j)))*sqrt(h)*Z1(j))...
                +(sigma*sqrt(1-rho^2)*sqrt(max(0, vec_V(j)))*sqrt(h)*Z2(j));
            vec_D(j+1) = vec_D(j)+((mu*vec_D(j))*h)+(vec_D(j)*sqrt(abs(vec_V(j)))*sqrt(h)*Z1(j));
            vec_X(j+1) = vec_X(j)+((r+(0.5*M_Var_con*vec_V(j))-(1+theta)*vec_X(j))*alpha*h)...
                +(alpha*M_sigma_1*sqrt(abs(vec_V(j)))*sqrt(h)*Z1(j))...
                +(alpha*M_sigma_2*sqrt(abs(vec_V(j)))*sqrt(h)*Z2(j));        
            
        end  %end for j
        
        vec_S= vec_D.*exp(M_a+(M_b*vec_V)+(M_c*vec_X)); %stock price
        path_ret=(vec_S(J)/vec_S(1))-1;       
        
        if path_ret>0.05       % Checking weather return is greater than 0.05
            St=real(vec_S(end));     % Current stock price
            vt=M_Var_con*vec_V(end);  % Current stock return variance is equal to its long run mean
            Xt=vec_X(end);
            Put_OTM = RN_OptionPrices(M_kappa_v,M_zeta_v,M_eta_v,M_sigma_v,M_rho_sv,Tt,K_OTM_P*St,r,St,vt,Xt,alpha,'P');
            Put_ATM = RN_OptionPrices(M_kappa_v,M_zeta_v,M_eta_v,M_sigma_v,M_rho_sv,Tt,K_ATM*St,r,St,vt,Xt,alpha,'P');
            Put_OTM_Pos_IV(path,th)= blsimpv(St, K_OTM_P*St, r, Tt, real(Put_OTM), 1, 0, [], {'Put'});
            Put_ATM_Pos_IV(path,th)= blsimpv(St, K_ATM*St, r, Tt, real(Put_ATM), 1, 0, [], {'Put'});
            
            Call_OTM = RN_OptionPrices(M_kappa_v,M_zeta_v,M_eta_v,M_sigma_v,M_rho_sv,Tt,K_OTM_C*St,r,St,vt,Xt,alpha,'C');
            Call_ATM = RN_OptionPrices(M_kappa_v,M_zeta_v,M_eta_v,M_sigma_v,M_rho_sv,Tt,K_ATM*St,r,St,vt,Xt,alpha,'C');
            Call_OTM_Pos_IV(path,th)= blsimpv(St, K_OTM_C*St, r, Tt, real(Call_OTM), 1, 0, [], {'Call'});
            Call_ATM_Pos_IV(path,th)= blsimpv(St, K_ATM*St, r, Tt, real(Call_ATM), 1, 0, [], {'Call'});
            
        end  %end for high positive return paths
        
        if path_ret<-0.05      % Checking weather return is less than 0.05
            St=real(vec_S(end)); % Current stock price
            vt=M_Var_con*vec_V(end);  % Current stock return variance is equal to its long run mean
            Xt=vec_X(end);
            Put_OTM = RN_OptionPrices(M_kappa_v,M_zeta_v,M_eta_v,M_sigma_v,M_rho_sv,Tt,K_OTM_P*St,r,St,vt,Xt,alpha,'P');
            Put_ATM = RN_OptionPrices(M_kappa_v,M_zeta_v,M_eta_v,M_sigma_v,M_rho_sv,Tt,K_ATM*St,r,St,vt,Xt,alpha,'P');
            Put_OTM_Neg_IV(path,th)= blsimpv(St, K_OTM_P*St, r, Tt, real(Put_OTM), 1, 0, [], {'Put'});
            Put_ATM_Neg_IV(path,th)= blsimpv(St, K_ATM*St, r, Tt, real(Put_ATM), 1, 0, [], {'Put'});
            
            Call_OTM = RN_OptionPrices(M_kappa_v,M_zeta_v,M_eta_v,M_sigma_v,M_rho_sv,Tt,K_OTM_C*St,r,St,vt,Xt,alpha,'C');
            Call_ATM = RN_OptionPrices(M_kappa_v,M_zeta_v,M_eta_v,M_sigma_v,M_rho_sv,Tt,K_ATM*St,r,St,vt,Xt,alpha,'C');
            Call_OTM_Neg_IV(path,th)= blsimpv(St, K_OTM_C*St, r, Tt, real(Call_OTM), 1, 0, [], {'Call'});
            Call_ATM_Neg_IV(path,th)= blsimpv(St, K_ATM*St, r, Tt, real(Call_ATM), 1, 0, [], {'Call'});
            
        end  %end for low negative return paths        
        
    end  % end for path simulation
    
end  % end for theta
toc
 

Put_OTM_Pos_mean=round(nanmean(Put_OTM_Pos_IV)*100,2)
Put_OTM_Neg_mean=round(nanmean(Put_OTM_Neg_IV)*100,2)
Put_ATM_Pos_mean=round(nanmean(Put_ATM_Pos_IV)*100,2)
Put_ATM_Neg_mean=round(nanmean(Put_ATM_Neg_IV)*100,2)

Call_OTM_Pos_mean=round(nanmean(Call_OTM_Pos_IV)*100,2)
Call_OTM_Neg_mean=round(nanmean(Call_OTM_Neg_IV)*100,2)
Call_ATM_Pos_mean=round(nanmean(Call_ATM_Pos_IV)*100,2)
Call_ATM_Neg_mean=round(nanmean(Call_ATM_Neg_IV)*100,2)
 
%% THE END

 
 
 
 
 
 
 